function a2i_measurements_archive, parameters, time, conc, sigma, pi, representation_error

  compile_opt idl2, hidden

  fname_search=a2i_filestr(/input, parameters['CASE_NAME'] + '/measurements/' + $
    parameters['POLLUTANTS', pi] + '_archive.csv')
  fname=file_search(fname_search)

  if fname ne '' then begin
  
    flines=file_lines(fname)-1
    archive_conc=fltarr(flines)
    archive_sigma=fltarr(flines)
    archive_box=intarr(flines)
    archive_time=fltarr(flines)
  
    openr, fun, fname, /get_lun
      str=''
      header=''
      n=0
      while not eof(fun) do begin
        readf, fun, str        
        if str ne '' and $
          strUpCase(strmid(strtrim(str, 2), 0, 4)) ne 'TIME' and $
          strmid(strtrim(str, 2), 0, 1) ne '%' and $
          strmid(strtrim(str, 2), 0, 2) ne '"%' and $
          strmid(strtrim(str, 2), 0, 1) ne ';' and $
          strmid(strtrim(str, 2), 0, 2) ne '";' then begin
          data=strSplit(str, ',', /extract)
          archive_time[n]=data[0]
          archive_box[n]=fix(data[1])
          archive_conc[n]=data[2]
          archive_sigma[n]=data[3]
          n++
        endif
      endwhile
    close, fun
    free_lun, fun
    
    wh=where(archive_time ne 0., nArchive)
    if nArchive gt 0 then begin
      archive_time=archive_time[wh]
      archive_box=archive_box[wh]
      archive_sigma=archive_sigma[wh]
      archive_conc=archive_conc[wh]
    endif
    
    caldat, time, m, d, y
    caldat, mr_timedec(archive_time), am, ad, ay

    for n=0, nArchive-1 do begin
      wh=where(m eq am[n] and y eq ay[n], count)
      if count gt 0 then begin
        conc[archive_box[n], wh]=mean([conc[archive_box[n], wh], archive_conc[n]], /nan)
        if finite(sigma[archive_box[n], wh]) then begin
          sigma[archive_box[n], wh]=sqrt(sigma[archive_box[n], wh]^2 + archive_sigma[n]^2)
        endif else begin

          wh_representation=where(finite(conc[archive_box[n], *]), count)
          if count eq 0 then begin
            message, 'REPRESENTATION ERROR IN ARCHIVE'
          endif

          sigma[archive_box[n], wh]=sqrt($
            2.*(mean(representation_error[archive_box[n], wh_representation[0]:-1]/conc[archive_box[n], wh_representation[0]:-1], /nan)*archive_conc[n])^2 + $
            archive_sigma[n]^2)
            
        endelse
      endif
    endfor
    
    wh=where(sigma eq 0., count)
    if count gt 0 then begin
      sigma[wh]=mean(sigma, /nan)
    endif
  
  endif

  out={conc:conc, sigma:sigma}
  return, out

end

; docformat = 'rst'
;
;+
;
; :Purpose:
;   Format AGAGE measurements for inversion
;
; :Inputs:
;   Parameters hash table
;   
; :Keywords:
;   random_scale_error: (boolean, input), set to 1 to add a random error to the measurement scale
;   
; :Outputs:
;   Measurements hash table
;
; :History:
; 	Written by: Matt Rigby, MIT, Aug 1, 2012
;
;-
function a2i_measurements, parameters, random_scale_error=random_scale_error

  compile_opt idl2, hidden

  y_pollutant=!null
  y_box=!null
  y_ti=!null
  y=!null
  y_error=!null
  y_start=intarr(parameters['NPOLLUTANTS'])
  y_first=fltarr(parameters['NPOLLUTANTS'])

  for pi=0, parameters['NPOLLUTANTS']-1 do begin

    precision_function=fltarr(parameters['TIMESIZE'])
    precision=float(strsplit(parameters['MEASUREMENT_REPEATABILITY', pi], $
      ';', /preserve_null, /extract))

    ;PRECISIONS
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    if n_elements(precision) gt 1 then begin
      precision_time=findgen(parameters['TIMESIZE'])/12. + parameters['STARTY']
      precision_year=float(strsplit(parameters['MEASUREMENT_REPEATABILITY_YEAR', pi], $
        ';', /preserve_null, /extract))
      for yi=0, n_elements(precision_year)-2 do begin
        wh=where(precision_time ge precision_year[yi] and precision_time lt precision_year[yi+1], count)
        if count gt 0 then begin
          precision_function[wh]=findgen(count)*(precision[yi+1] - precision[yi])/float(count) + $
            precision[yi]
        endif
      endfor
    endif
    
    if n_elements(precision) eq 1 then begin
      precision_function[*]=precision[0]
    endif
    
    precision_function=precision_function/100.

    ;Get AGAGE DATA
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ;DATES TO REMOVE
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    remove=strsplit(parameters['MEASUREMENT_REMOVE', pi], ';', /preserve_null, /extract)
    
    exclude=!null    

    if remove[0] ne '' then begin

      exclude=hash()

      for ri=0, n_elements(remove)-1 do begin

        remove_station=strTrim(remove[ri], 2)
        key=strUpcase(strmid(remove_station, 0, 3))
        remove_station=strTrim(strMid(remove_station, 3), 2)
        remove_station=strMid(remove_station, 1, strLen(remove_station)-2)

        remove_periods=strSplit(remove_station, ':', /extract)
        
        exclude_station=dblarr(n_elements(remove_periods)*2)

        for period=0, n_elements(remove_periods)-1 do begin
          remove_period=strSplit(remove_periods[period], ' to ', /extract)
          if n_elements(remove_period) ne 2 then message, "'MEASUREMENTS TO REMOVE' MUST INCLUDE 'TO'"
          exclude_station[period*2]=time_convert(remove_period[0], 'M/D/Y')
          exclude_station[period*2+1]=time_convert(remove_period[1], 'M/D/Y')
        endfor
        
        exclude[key]=exclude_station

      endfor

    endif
    
    
    ;GET MEASUREMENTS
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    agage_get_box, time, conc, pollutant=parameters['POLLUTANTS', pi], $
      startT=julday(1, 1, parameters['STARTY'], 0), EndT=julday(1, 1, parameters['ENDY'], 0), $
      error_variability=variability_AGAGE, error_precision=precision_agage, $
      exclude=exclude, ale=parameters['ALE', pi]
    sigma=fltarr(4, parameters['TIMESIZE'])

    if n_elements(conc) eq 0 then begin
      Alternative_measurements=1
    endif else begin
      Alternative_measurements=0
    endelse

    ;REPLACE WITH MEASUREMENTS FROM FILE, IF REQUIRED
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    if parameters['MEASUREMENT_FOLDER'] ne '' or Alternative_measurements then begin
      if Alternative_measurements then begin
        fname=file_search(a2i_filestr(/Input, parameters['CASE_NAME'] + '/measurements_alt/' + $
          parameters['POLLUTANTS', pi] + '.sav'))
        if fname eq '' then begin
          message, 'Put some ' + parameters["POLLUTANTS", pi] + ' measurements in measurements_alt/'
        endif
        sObj=OBJ_new('IDL_Savefile',fname)
        sVarNames=sObj->names()
        restore, fname
      endif else begin
        sObj=OBJ_new('IDL_Savefile', a2i_filestr(/Input, parameters['CASE_NAME'] + '/' + $
          parameters['MEASUREMENT_FOLDER'] + '/' + parameters['POLLUTANTS', pi] + '.sav'))
        sVarNames=sObj->names()
        for vi=0, n_elements(sVarNames)-1 do sObj->restore, sVarNames[vi]
      endelse
      
      wh=where(time ge julday(1, 1, parameters['STARTY'], 0) and time lt julday(1, 1, parameters['ENDY'], 0), count)
      if count eq 0 then begin
        message, strcompress('No measurement file: ' + parameters['POLLUTANTS', pi] + '.sav')
      endif else begin
        time_temp=timegen(parameters['TIMESIZE'], start=julday(1, 15, parameters['STARTY'], 0), units='MONTHS')
        caldat, time_temp, month_temp, day_temp, year_temp
        caldat, time, month, day, year
        conc_temp=replicate(!values.f_nan, 4, parameters['TIMESIZE'])
        sigma_temp=replicate(!values.f_nan, 4, parameters['TIMESIZE'])
        for ti=0, parameters['TIMESIZE']-1 do begin
          wh=where(month eq month_temp[ti] and year eq year_temp[ti], count)
          if count gt 0 then begin
            conc_temp[*, ti]=conc[*, wh[0]]
            ;if uncertainty given in file, overwrite existing sigma
            if total(sVarNames eq 'SIGMA') gt 0 then begin
              sigma_temp[*, ti]=sigma[*, wh[0]]
            endif
          endif
        endfor
      endelse

      conc=temporary(conc_temp)
      sigma=temporary(sigma_temp)
      time=temporary(time_temp)

    endif

    if strmatch(parameters['CASE_NAME'], 'NOAA') and total(sigma) gt 0. then begin
      print, "WARNING, MAKE SURE MEAUSREMENT UNCERTAINTY ISN'T IN MEASURMENT FILE"
    endif

    ;CALCULATE UNCERTAINTY
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    if total(sigma) eq 0. or total(finite(sigma)) eq 0 then begin
      for bi=0, 3 do begin
        sigma[bi, *]=sqrt(variability_AGAGE[bi, *]^2 + (precision_function*conc[bi, *])^2)
        wh=where(finite(conc[bi, *]) and (finite(variability_AGAGE[bi, *]) eq 0), count)
        if count gt 0 then sigma[bi, wh]=sqrt(mean(variability_AGAGE[bi, *], /nan)^2 + $
          (precision_function*conc[bi, wh])^2)
      endfor
    endif

    ;Include additional measurements (optional)
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    if parameters.haskey('ARCHIVE') then begin
    if parameters['ARCHIVE'] then begin
      out=a2i_measurements_archive(parameters, time, conc, sigma, pi, variability_AGAGE)
      conc=out.conc
      sigma=out.sigma
    endif
    endif

    ;Include detection limit (optional)
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    whMin=where(sigma lt parameters['MEASUREMENT_DETECTION_LIMIT', pi], count)
    if count gt 0 then begin
      sigma[whMin]=parameters['MEASUREMENT_DETECTION_LIMIT', pi]
    endif

    ;Apply R1 scaling (optional)
    if keyword_set(random_scale_error) then begin
      conc=(1. + randomn(systime(/seconds)))*parameters['SCALE_ERROR_PLUS', pi]*conc
    endif

    ;PREPARE OUTPUT VECTORS
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    caldat, time, dummy, dummy1, year
    wh=where(year lt parameters['POLLUTANT_START_YEAR', pi], count)
    if count gt 0 then begin
      conc[*, wh]=!values.f_nan
    endif
    
    
    y_exists=intarr(n_elements(conc[0, *]))
    for bi=0, 3 do y_exists+=finite(conc[bi, *])

    wh_exists=where(y_exists gt 0, count_exists)
    
    if ~parameters['FIRN'] then begin
      y_first[pi]=mean(conc[*, wh_exists[0]], /nan)
    endif

    if count_exists gt 0 then begin
      if ~parameters['FIRN'] then begin
        y_start[pi]=12*(wh_exists[0]/12)
      endif

      for bi=0, 3 do begin
        y=[y, reform(conc[bi, y_start[pi]:-1])]
        y_error=[y_error, reform(sigma[bi, y_start[pi]:-1])]
        y_box=[y_box, replicate(bi, n_elements(time) - y_start[pi])]
        y_ti=[y_ti, indgen(n_elements(time) - y_start[pi]) + y_start[pi]]
      endfor
      y_pollutant=[y_pollutant, replicate(pi, (n_elements(time) - y_start[pi])*4)]
    endif else begin
      print, 'Warning, no measurements: ', parameters['POLLUTANTS', pi]
    endelse

  endfor  ;pollutant


  ;IDENTIFY MISSING MEASUREMENTS
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  wh_zero=where(y le 0., count)
  if count gt 0 then begin
    y[wh_zero]=!values.f_nan
  endif
  
  wh_measure=where(finite(y), count)
  if count eq 0 then begin
    message, 'NO MEASUREMENTS'
  endif


  ;Apply function
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  for pi=0, parameters['NPOLLUTANTS']-1 do begin
    wh=where(y_pollutant eq pi)
    case strUpCase(parameters['FUNCTION_Y', pi]) of
      'LN(Y)': begin
        y_error[wh]=y_error[wh]/y[wh]
        y[wh]=alog(y[wh])
      end
      'Y': ;Do nothing
    endcase
  endfor

  ;OUTPUT
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  output=hash()
  output['Y']=y
  output['POLLUTANT']=y_pollutant
  output['BOX']=y_box
  output['TI']=y_ti
  output['ERROR']=y_error
  output['START']=y_start
  output['FIRST']=y_first
  output['NMEASUREMENTS']=n_elements(y)
  output['WH_MEASURE']=wh_measure

  return, output

end